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The observed angular power spectrum of the cosmic microwave background temperature 
anisotropy, Cj, is a convolution of a cosmological radiative transport kernel with an assumed pri- 
mordial power spectrum of inhomogeneities. Exquisite measurements of Ci over a wide range of 
multipoles from the Wilkinson Microwave Anisotropy Probe (WMAP) has opened up the possibil- 
ity to deconvolve the primordial power spectrum for a given set of cosmological parameters (base 
model). We implement an improved (error sensitive) Richardson-Lucy deconvolution algorithm on 
the measured angular power spectrum from WMAP assuming the concordance cosmological model. 
The most prominent feature of the recovered P{k) is a sharp, infra-red cut off on the horizon scale. 
The resultant Ci spectrum using the recovered spectrum has a likelihood far better than a scale 
invariant, or, 'best fit' scale free spectra {AlnC ~ 25 w.r.t. Harrison Zeldovich, and, Aln£ ~ 11 
w.r.t. power law with Us = 0.95). The recovered P{k) has a localized excess just above the cut-ofi^ 
which leads to great improvement of likelihood over the simple monotonic forms of model infra-red 
cut-off spectra considered in the post WMAP literature. The recovered P{k), in particular, the 
form of infra-red cut-off is robust to small changes in the cosmological parameters. We show that 
remarkably similar form of infra-red cutoff is known to arise in very reasonable extensions and re- 
finements of the predictions from simple infiationary scenarios. Our method can be extended to 
other cosmological observations such as the measured matter power spectrum and, in particular, 
the much awaited polarization spectrum from WMAP. 

PACS numbers: 98.80.Es, 98.70.Vc 



I. INTRODUCTION 



Increasingly accurate measurements of the anisotropy in the temperature of the cosmic mi- 
crowave background (CMB) has ushered in an era of precision cosmology. A golden decade of 
CMB anisotropy measurements by numerous experiments was topped by the results from the first 
year of data obtained by the Wilkinson Microwave Anisotropy Probe (WMAP) 0|. Under simple 
hypotheses for the spectrum of primordial perturbations, exquisite estimates of the cosmologi- 
cal parameters have been obtained from the angular power spectrum measurement by WMAP 
combined with other cosmological observations Although the assumed, scale free (with mild 
deviations) , initial power spectra may be a generic prediction of the simplest scenarios of gener- 
ation of perturbations during inflation, initial spectra with radical deviations are known to arise 
from very reasonable extensions, or, reflnements to the simplest scenarios Consequently, 
cosmological parameter estimation from the CMB anisotropy and the matter power spectrum ob- 
tained from redshift surveys, weak gravitational lensing and Ly-a absorption, depends sensitively 
on the dimensionality, nature and freedom in the parameter space of initial conditions 

The angular power spectrum, C;, is a convolution of the initial power spectrum P{k) generated 
in the early universe with a radiative transport kernel, G{1, k), that is determined by the current 
values of the cosmological parameters. The remarkably precise observations of the angular power 
spectrum Ci by WMAP, and the concordance of cosmological parameters measured from different 
cosmological observations opens up the avenue to directly recover the initial power spectrum of 
the density perturbation from the observations. The Richardson-Lucy (RL) method deconvolution 
was shown to be a promising and powerful method to measure the power spectrum of initial 
perturbations from the CMB angular power spectrum 01 ■ In this paper, we apply the method 
to the CMB anisotropy spectrum measured by WMAP. We have also devised and implemented 
an improvement to the RL scheme, whereby the iterative deconvolution algorithm is designed to 
converge and match the measurements only within the given error-bars. 

The direct numerical deconvolution has clear advantages over the other prevalent approach of 
obtaining the most likely parameter values for a parametric model primordial spectra [52]. First, 
as emerges from our work, the direct method can reveal features that are not anticipated by the 
parametric model spectra, hence would be completely missed out in the latter approach. Second, in 
the absence of an accepted early universe scenario (more narrowly, a favored model of inflation) , it is 
difficult to a priori set up and justify the chosen space of initial conditions. The complex covariances 
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between the cosmological and the initial parameters are sensitive to the parameterization of the 
space of initial spectra adopted. Efforts along these lines are further obscured by issues such 
as the applicability of the Occam's razor to dissuade extension of the parameter space of initial 
conditions. Such deliberations have been recently framed in the more quantitative language of 
Bayesian evidence to evaluate and select between possible parameterizations However, this 

approach cannot really point to a preferred parameterization. Whereas, in the direct approach we 
can evade the issue of appropriate parameterization of the initial power spectrum. For a given set 
of cosmological parameters, our method obtains the primordial power spectrum that 'maximizes' 
the likelihood. Hence, the space of parameters remains that of the (more widely accepted and 
agreed up on) cosmological parameters alone. In principle, it is possible to explore the entire space 
of cosmological parameters along the lines being done routinely, and instead of simply computing 
the likelihood for a given model initial power spectrum, one obtains the initial power spectrum 
that maximizes the likelihood at that point and assigns that likelihood to that point in the space of 
cosmological parameters However, in this paper we limit ourselves to obtaining the primordial 
spectrum for a few given sets of 'best fit' cosmological parameters and showing that recovered P{k) 
is robust to small local variations in the cosmological parameters. We defer the wide exploration 
of the cosmological parameter space to future work. 

In section ^ we describe the problem and present the deconvolution method, details of its 
implementation and tests of recovering known primordial spectra from synthetic angular power 
spectrum data. The recovered spectrum from WMAP data using our method is described in 
section UTTl The recovered primordial spectrum has statistically significant features that are robust 
to the variations in the cosmological parameters. In section Hvl the recovered spectrum is shown 
to be tantalizingly similar to some forms of broken scale invariant spectra known to arise from 
reasonable variations of the simplest inflationary scenario. 

II. METHOD 

A. Integral equation for CMB anisotropy 

In this subsection we recall the integral equation for the angular power spectrum of CMB 
anisotropy and set up the inverse problem that we solve using a deconvolution method. The 
observed CMB anisotropy AT(n) is one realization of a random field on the surface of a sphere 
and can be expressed in terms of the random variates aim that are the coefficients of a Spherical 
Harmonic expansion given by 

oo m—l 

AT(n) = X! ai,nYim{n) (2.1) 

1 — 2 ni— — l 

where y/m(n) are the Spherical Harmonic functions. For an underlying isotropic, Gaussian statis- 
tics, the angular power spectrum, Ci defined through 

{airnavm') = ClSll'S„im' , (2.2) 

completely characterizes the CMB anisotropy. 

In a flat universe the temperature fluctuation in the CMB photons at the location x at the 
present conformal time 770 propagating in a direction n is 

A(n)= A(xo,n,r;o) = J d^k e^'^-^ A{k,n,m). (2.3) 

For globally isotropic cosmology, the temperature fluctuation A(k, 11,770) = A(k • 11,770) can be 
expanded in terms of Legendre polynomials leading to 



AT(2:,n,77o) = / rf^A^e^"-^ ;^(-z)'(2; + l)Ai(fc,77o)fl(k.ii). (2.4) 

1=0 
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The angular power spectrum Ci given by the coefficients of Legendre expansion is then expressed 

as 



/rik 
^P(fc)|AT/(fc,r7o)P (2.5) 

where P{k), the power spectrum of the primordial (scalar) metric initial perturbation V'prim, is 
given by 



(V'pnm(k) V;rim(k')) ^ ^ 6{k - k') (2.6) 

because the k space modes are uncorrelated in a homogeneous space. The spectrum P{k) represents 
the r.m.s. power in the scalar metric perturbations per logarithmic interval dk/k at a wavenumber 
k. It is related to power spectrum of the primordial modes of density perturbations, Sk, as P{k) cx 
|(5fep//c. For the conventional scale free parameterization, the spectral index Ug is defined through 
|(5fcp — Ak"^" . The power spectrum P{k) is a constant for the scale invariant Harrison- Zeldovich 
spectrum (corresponds to ris = 1). 

The harmonics of the temperature fluctuations at the current epoch, ATi{k, r]o), is obtained from 
the solution to the Boltzmann equation for the CMB photon distribution. In this paper we use 
ATi{k,r]o) computed by the CMBfast software Numerically, a suitably discretized space of 

wave-numbers, ki is used where the following discrete version of integral eq. Il2.5|l is applicable 



Ci =^G(;,A;,)P(fe^) 

i 

Ak- 

G{l,h) = —^\ATi{h,7jo)\^. (2.7) 

In the above equation, the 'target' angular power spectrum, C; = Cf , is the data given by obser- 
vations, and the radiative transport kernel, G{1, k) is fixed by the cosmological parameters of the 
'base' model. (The kernel G{1, k) also includes the effect of geometrical projection from the three 
dimensional wavenumber, /c, to the harmonic multipole, I on the two dimensional sphere.) For a 
given G{l,k), obtaining the primordial power spectrum, P{k) from the measured C/ is clearly a 
deconvolution problem. An important feature of our problem is that Cf , G{1, k) and P{k) are all 
positive definite. However, to get reliable results from the deconvolution, we require high signal 
to noise measurements of over a large range of multipoles with good resolution in multipole, 
preferably, from a single experiment to avoid the uncertainties of relative caHbration 0| . 

Ideally, it will be best to measure each Gf independently. In practice, incomplete sky coverage 
and other effects limit the resolution in multipole space. All experiments provide band power 
estimates, Cb{lcii) = J^i^^i^^ which are averaged Hnear combinations of the underlying Gi. 
Since eq. 112.711 is linear, a similar equation holds for all the band powers with a kernel, G{lcs, h) = 

w/^"* G(Z, fci). For simplicity and brevity of notation, we retain the notation Gf and G{l,k) 
for band power estimates with I denoting the bin center. We also implicitly assume a discrete 
wavenumber k instead of carrying the clumsy notation, ki . 

As mentioned in more detail in section UlI AL the CMB angular power spectrum from WMAP 
ranges from the quadrupole, Z = 2 to 900. Moreover, the 'full sky' coverage of WMAP impHes 
good resolution in multipole space. We use WMAP TT (temperature-temperature) binned power 
spectrum as Gf in eq. H2.7|l . We use CMBfast software to compute the the G{1, k) matrix for the 
post- WMAP 'best fit' cosmological parameters. We emphasize that although the P{k) is recovered 
using binned data, the significance (performance) of the recovered spectrum is evaluated using 
the WMAP likelihood of entire unbinned C/ properly accounting for covariances. The likelihood 
is computed using the numerical code, data and its error covariance provided by the WMAP 
collaboration with the release of the first year data. 

In this work we limit our attention to the angular spectrum of the temperature anisotropy, Cf^. 
Including the polarization of CMB photons, equations similar to eq. 112. 5|l can be also written for 
the three additional angular power spectra, Cf^, Gf^ and G^^ involving their corresponding 
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kernels. It will certainly be interesting to include these when more complete polarization data is 
available in the future. At present, only the Cf^ spectrum has been published by the WMAP 
team which is not positive definite and hence, not ideally suited for our method. However, once 
Cf'^ data is made available, our method can readily accommodate both Cf^ and Cf"^ together 
since combinations such as Cf^ ± 2Cj^^ + Cj^^ are positive definite. 



B. Deconvolution method 

The Richardson-Lucy (RL) alg orithm was developed and is widely used in the context of image 
reconstruction in astronomy fl^ . However, the method has also been successfully used in 
cosmology, to deproject the 3-D correlation function and power spectrum from the measured 2-D 
angular correlation and 2-D power spectrum [18, 19]. We employ an improved RL method to solve 
the inverse problem for P{k) in eq. 112.711 . The advantage of RL method is that positivity of the 
recovered P{k) is automatically ensured, given G{l,k) is positive definite and C;'s are positive. 

The RL method is readily derived from elementary probability theory on distributions To 
make this connection, we consider normalized quantities 

J2Ci^l; J2 ^(^) ^ 1; E ^('' ^) = 1 ■ (2.8) 

/ k I 

This allows us to view the functions, P{k) and Ci as one dimensional distributions, and G{1, k) 
as a conditional probability distribution. Further, for convenience (not necessity) of writing an 
infinitesimal measure dl we view I to be continuous. The integrand in the integral eq. 112. 5|l suggests 
defining two other probability distributions Q{l,k) and L{l,k), such that 

Gil, k) P{k) dldk = Q{1, k) dldk = L{k, I) Ci dldk (2.9) 
Dividing the both side of eq. (I2.9|l by C; dldk we obtain 

mi) = ^G{l,k). (2.10) 

The normalization conditions imply 

P{k) = J Qil,k)dl = J CiL{k,l)dl, (2.11) 
which in the discrete I space reads 

P(fc) = ^) = E ^) ■ (2.12) 

The RL method iteratively solves the eqs. II2.1()|I and H2.12II . Starting from an initial guess 
P(°)(fc), L is obtained using eq. Ij2.1fl|l as the first step. The second step is to obtain a revised 
p(^^(fc) using eq. II2.12|I . These two steps are repeated iteratively with the P'^^\k) obtained after 
iteration i feeding into the iteration i + \. In principle, the final answer could depend on the 
initial guess but in practice, for a large variety of problems, RL is known leads to the correct 
answer even a crude estimation of the initial guess. In particular, for our problem the RL rapidly 
converges to the same solution P{k) independent of the initial guess p(°)(fc). This is demonstrated 
in Annendix IA II 

The iterative method can be neatly encoded into a simple recurrence relation. The power spec- 
trum p(*+^)(fc) recovered after iteration (i -I- 1) is given by 

fD _ Mi) 

P^'+^Hk) - P^'Hk) = pW(fc) J2 G{l,k)^ (-p^- (2.13) 

I 
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where Cf is the measured data (target) and C; is the angular power spectrum at z*'' iteration 
obtained from eq. I|2.7|l using the recovered power spectrum P^^\k). Eq. H2.13II with eq. 112.711 for 
obtaining C^^^ from P*-*-* (fc) completely summarizes the standard RL method. 

Due to noise and sample variance, the data Cf is measured within some non-zero error bars ai. 
The standard RL method does not incorporate the error information at all. Consequently, a well 
known drawback of the standard RL method is that at large iterations the method starts fitting 
features from the noise. Modified forms of RL that address this issue have been proposed (e.g., 
see damped RL method in [ISl)- In our problem, this problem manifests itself as very non-smooth 
deconvolved spectrum P{k) that has poor likeHhood with the full WMAP spectrum data. We 
devise a novel method to make the RL method sensitive to the error ai by modifying eq. H2.13II to 



p^'+^Hk) - p^'Hk) = p«(fc) y Gil, k) 



Ci -c. 



{i) 



I 



tanh 



(2.14) 



The idea is to employ a 'convergence' function to progressively weigh down the contribution to the 
correction p(*+i) — p(') from a multipole bin where c[^^ is close Cj^ within the error bar ct;. This 
innovation significantly improves the WMAP likelihood of the deconvolved spectrum. For certain 
G(l, k), the improvement is so dramatic that using IRL becomes very crucial to successful recovery 
of the spectrum (see IIII Ct . The final results are not sensitive to the exact functional form of the 
convergence function. The choice given in eq. H2.14II works well but is not unique in any sense. 

At every iteration of the IRL scheme, we compute the of the cj;^^ with respect to the binned 
data Cf . We have found that the IRL iterations (as well as the RL) march almost monotonically 
toward improved (smaller) ■ We halt the iterations when the does not change appreciably in 
subsequent iterations of IRL. 
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FIG. 1: The curves are G{l,k) versus wavenumber k used in our work. G{l,k) is averaged over G{l,k) 
within multipole bins used by WMAP. The two vertical lines roughly enclose the region of fc-space strongly 
probed by the kernel where the primordial spectrum can be expected to be recovered reliably. The fc-space 
sampling used is indicated by the line of ' + ' symbols at the top of the plot. 



C. Post processing the deconvolved spectrum 



The deconvolution algorithm produces a 'raw' P(fc) that has to be processed further. The raw 
deconvolved spectrum has spurious oscillations and features arising largely out of the k space 
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FIG. 2: Each of the three rows in the panel of figures illustrates the recovery the primordial power for a test 
case using Ci arising from a known non-scale invariant primordial spectrum. The first column compares the 
raw deconvolved spectrum with the input spectrum. Note the similar artifacts in the all the raw spectra at 
the low and high k end discussed in the text. The feature is outside the range of G(i, k) and is completely 
missed in the third case. The second, column is the differenced spectrum obtained by dividing out by 
a raw reference spectrum. The differenced spectrum resembles the input spectrum (top two cases) with 
small oscillations. The third column shows that the final recovered spectrum obtained by smoothing the 
differenced spectrum matched the input spectrum very well. 



sampling and binning in I space. These numerical artifacts are common to all recovered spectra 
recovered with same Gil, k). Fig.n]shows the plot of the binned kernel G(/, k) for the /-space bins 
used by WMAP ES). The wavenumbers are scaled by horizon scale, kh — 27r/(?7o — jyrcc)- The 
sampling of k space used here is also indicated in the figure. We find that removing these artifacts 
and smoothing the resultant spectrum improves the WMAP likeHhood of the corresponding Ci . 

We illustrate the steps of removing the numerical artefacts and smoothing with test case exam- 
ples of synthetic C/ from known test primordial spectrum. The first column of Fig. 12 shows the 
'raw' spectrum obtained. A comparison with the known spectra shows similar numerical noise and 
artefacts in all the cases, especially, the rise at very low and very high wave numbers. As shown 
in a.nnendix lA 31 the main artifacts at the low and high k ends can be understood and modeled 
analytically. However we find it easier to remove them constructing a numerically generated tem- 
plate reference spectrum that also takes care of small features which appear due to changes in 
the k spacing. To remove these numerical artifacts, we generate a synthetic C/ using a Harrison- 
Zeldovich spectrum and then apply the deconvolution method to recover a 'reference' spectrum 
Pref(^)- We divide every 'raw' spectrum by Pref(fc) to obtain a 'differenced' spectrum that does 
not have numerical artifacts as seen in middle column ("diff.") of Fig|2| 
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The difference spectra are a noisy version of the test spectra. Hence, the exercise with the test 
spectra suggests that the differenced spectrum needs to be 'suitably' smoothed to recover the true 
power spectrum. A simple smoothing procedure with simple fixed width window function leads to 
a satisfactory recovery of the test spectrum. (A 'Bowler hat' window constructed with hyperbolic 
tangent functions is used for smoothing.) 

The last column in Fig |2l shows the remarkably successful recovery of the test spectra in the top 
two case. As expected, the P{k) is recovered best in the k space (within the vertical Hues in bottom 
row of FigEJ where the kernel G{l,k) is significant (marked by the vertical lines in Fig.n]and[21l. 
In the bottom row, the feature in the test spectrum is outside the range of the kernel G{1, k), and, 
as expected, the recovery process misses it completely. For the top row, if we ignore the region 
where there is no power for G{1, k) the recovered spectrum is well matched with the test spectrum. 
The match may be further improved by using more elaborate, adaptive smoothing procedure. 

For application to real data, the phrase 'suitable smoothing procedure' may appear ambiguous. 
But the smoothing is in fact very well defined for real data by demanding that the smoothed P{k) 
produces a theoretical C; that has higher Hkelihood given the data. We find that this approach 
works extremely well. The deconvolution algorithm uses the binned C; data, hence the theoretical 
Ci corresponding to the recovered P{k) fits the binned data very well. However, the WMAP 
Hkelihood of the theoretical Ci suffers owing to spurious oscillations in the differenced spectra. 
The WMAP likelihood improves as the differenced spectra, Pdiff(fc) is smoothed. We smooth the 
differenced spectrum so that WMAP likelihood of the corresponding theoretical C; is maximized. 
Although it is difficult to establish that the final result is the unique solution with maximum 
Hkelihood, in practice, our simple scheme does lead to a well defined result (no distinct degenerate 
solutions were found). Since our smoothing procedure is simple minded, possible avenues for 
improvement with more elaborate smoothing procedure remain open. Work is in progress to 
employ wavelet decomposition for the smoothing procedure. 

III. APPLICATION TO THE WMAP CMB ANISOTROPY SPECTRUM 

We apply the method described in the previous section to the angular power spectrum obtained 
with the first year of WMAP data pubHcly released in February 2003 0| to recover the primordial 
power spectrum. In section lTlI Al the publicly available WMAP data and how it is used in our work 
is discussed. We also describe the choice of the 'base' cosmological model. The recovered primordial 
power spectrum for this model is presented in section ITlI Bl The next section ITlI CI presents the 
effect of varying the cosmological parameters (within 1 a error bars) on the recovered primordial 
power spectrum. We also present the primordial spectrum for a set of cosmological models with 
large optical depths (r = 0.1,0.17,0.25) corresponding to possible the early reionization scenarios 
suggested by the WMAP temperature-polarization (TE) cross-spectrum. 

A. WMAP anisotropy data and the cosmological model 

Accurate measurements of the angular power spectrum of CMB anisotropy was derived from the 
first year WMAP data recently [2l|. The spectrum obtained by averaging over 28 cross-channel 
power measurements is essentially independent of the noise properties of individual radiometers. 
The power at each multipole ranging from ^ = 20 to 900 was estimated together with the covari- 
ance The instrumental errors are smaHer than the cosmic variance up to Z ~ 350, and the 

signal to noise per mode is above unity up to ? ~ 650. 

The angular power spectrum estimate and covariance matrix are publicly available at the 
LAMBDA data archive The WMAP team has also made available a suite of F90 codes 

that computes the likelihood for a given theoretical C/ spectrum given the full angular power spec- 
trum included the covariance measured by WMAP. We use the TT likelihood code for computing 
the Hkelihood of Ci obtained from the recovered power spectrum and refer to these numbers as the 
'WMAP likeHhood' in the paper. 

In addition, the WMAP team has also obtained a binned angular power spectrum where an 
average Ci is defined over bins in multipole space. The binned C; estimates can be treated as 
independent data points since the covariance between binned estimates is negHgible. We use this 
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binned spectrum as in the deconvolution of eq. 112. 7|l . (We are aware of but do not consider the 
revised estimates of the low multipoles made by other authors after the WMAP results ^23., .24]) 
The variance of C/ measurements is given by 

where crjv is the noise per pixel, Wp the angular pixel size, /sky is the fraction of sky covered and 
Bi is the transform of the experimental beam [23|. It is important to note that contribution to 
the error from cosmic variance is proportional to the underlying theoretical/true Cf spectrum. 
To obtain the total error bars ai used in the IRL deconvolution eq. H2.14II . we should add the 
cosmic variance for the theoretical C; to the the statistical error bars given with the binned data. 
However, c[^^ rapidly iterates to within the error bars in the IRL method and the simpler 

option of using instead c[^^ for computing the cosmic variance works equally well since ai in 
eq. II2.14|I simply regulates the convergence of IRL [5^ . 

We consider a flat A-CDM universe, with Hubble constant, = 71 km.s~^/Mpc., Baryon 
density, 0.^ = 0.0224, a cosmological constant corresponding to f^A = 0.73, with the remaining 
balance of matter to critical density in cold dark matter. This is the 'concordance' cosmological 
model suggested by the WMAP parameter estimation. While we mainly focus attention on a 
cosmological model without early reionization (optical depth to reionization, t = 0) we do present 
in section UlI CI the recovered primordial spectra for models with early reionization with opacity 
going up to T = 0.25. The case for a large optical depth t comes from the temperature-polarization 
cross correlation. For simplicity, we Hmit ourselves to the temperature anisotropy spectrum and 
avoid undue attention to the cosmology suggested by the yet incomplete polarization data. The 
polarization spectrum is expected to be announced by WMAP soon. It is then easy to extend our 
method to include both the temperature and polarization anisotropy spectra. 
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FIG. 3: The three stages leading to the final recovered spectrum for a base cosmological model (r = 0.0, 
h = 0.71, nth^ = 0.0224 and = 0.73) is shown. The lower dashed line is the raw deconvolved power. 
The upper dashed line is the differenced spectrum obtained by dividing out by the reference spectrum 
shown as a dotted line. The solid line is the final result after smoothing that gives the best likelihood. 
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FIG. 4: The recovered Ci corresponding to the raw P{k) are shown in the upper row and that corresponding 
to the final smoothed P{k) spectrum are shown in the lower row. The left panels show the full range of 
multipole, while the right hand panels zoom into the low multipoles. The Cj's from the raw as well as 
final P{k) fit the binned Cf well (x^ ~ 10 and X ~ 20, respectively for 38 points). However, the jagged 
form of the C; between the I bins (apparent in the low multipoles on the right) leads to poor WMAP 
likelihood for the Ci from the raw P{k). The differencing and smoothing procedure irons out the jagged 
Ci dramatically improving the WMAP likelihood for the final smoothed P{k). The WMAP likelihood is 
more relevant since it incorporates the estimation of each Ci and the full error covariance. 



B. Primordial power spectrum from WMAP 

We apply the method described in the previous section to the WMAP data. Fig-Elshows the raw 
deconvolved spectrum, the reference spectrum, the differenced spectrum and the final recovered 
spectrum after smoothing obtained at each step of our method outlined in section Hll 

We discuss the advantage of using improved Richardson-Lucy (IRL) in appendix I A 21 The effect 
of IRL method is also evident in Fig. 31 The Ci from the recovered P{k) matches the binned CP 
only within the error-bars. The Ci spectrum corresponding to the differenced and smoothed P{k) 
is shown in Fig. 31 Although, the former has better for the binned data, the WMAP likelihood 
of the latter is much better. The poor likelihood of the differenced spectrum arises because there 
is no check on fluctuations in Ci at intermediate multipoles between the bin centers arising from 
the spurious numerical effects of /c-space and ^-space sampHng. Better WMAP likeHhood is more 
relevant than good for the binned data since the former incorporates the estimation of each 
Ci and the full error covariance. Thus the smoothing step our method is carried out with a well 
defined goal of maximizing the WMAP likelihood. 

The primordial spectra recovered from WMAP data is again shown in Fig.|3 The dark solid line 
is the primordial spectrum that has the best WMAP log-likelihood of In £ = —478.20 for the base 
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FIG. 5: The final recovered spectrum for the base cosmological model (r = 0.0, h = 0.71, = 0.0224 

and JIa = 0.73) is compared with set of P{k) with WMAP likelihood within ~ 2 u. The thick line gives 
the best likelihood equal to —478.2 and the other lines gives the likelihood bigger than —480. We can see 
that the sharp infra-red cut off is common to all these recovered spectra. The infra-red cutoff is remarkably 
close to the horizon scale and appears to be a robust feature. Another significant and robust feature is the 
bump just above the cut-off (reminiscent of the oscillation from under-damped transient). The difference 
between these spectra are in the smoothing and removing the noises from the raw deconvolved spectrum. 



cosmological model described in the previous section. The other lines have Hkelihood In £ > —480, 
i.e., roughly within 2 cr of the best one. For comparison, the same cosmological model with a 
scale invariant (Harrison-Zeldovich) primordial spectrum has In £ = —503.6, and, with a tilted 
scale free primordial spectrum n = 0.95, the likelihood improves to In £ = —489.3. (A comparison 
of likelihood numbers for various primordial spectra is given in Table The improvement in 
Hkelihood for the recovered spectra is striking. 

The most prominent feature of the recovered spectrum is the infra-red cutoff in the power 
spectrum remarkably close to the horizon scale (fc^ ~ fc/i = 27r/ryo) chiefly in response to the low 
quadrupole measured by WMAP. Another notable point is the slight tilt (n ~ 0.95) of the plateau 
in the recovered spectrum at large k which is consistent with the best fit n obtained by WMAP 
for power law primordial spectra. After the WMAP results, model power spectra with infra-red 
cutoffs of the form 



P{k) = A,k^ 



(3.2) 



invoked to explain the suppressed low multipoles of WMAP have had Hmited success (23, l27l |. 
While the effective xls = — 21n£ improved by at least 22 over a power law model, the model infra- 
red cutoff eq. Ij3.2ll could improve the XcS merely by ~ 3 [23|. Figure El compares the recovered 
spectrum with these model infra-red spectra. The power law spectrum shown in the figure also 
highlights the small tilt (n ^ 0.95) recovered by our method. 



C. Dependence on the cosmological parameters 



The primordial spectrum is recovered for a set of 'best fit' cosmological parameters defining the 
'base' model. We have varied each cosmological parameter by the quoted 1 a error-bars of the 
WMAP estimates [2j. Figure shows the dependence of the recovered power spectrum on small 
changes to cosmological parameters varied one at a time keeping the others fixed at their central 
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FIG. 6: Comparison between our recovered P(k) and exponential form of infra-red cut off recently studied 
to explain the suppressed quadrupole of WMAP angular power spectrum |2^,[23|- We note that the infra- 
red cut-off of the recovered spectrum is very steep. The power excess just above the cutoff in the recovered 
is extremely significant in the remarkably enhanced likelihood. The power spectrum (n = 0.95) shows that 
our method does recover the preferred tilt obtained by WMAP team parameter estimation with power law 
spectrum Q. 



value. We find that recovered spectrum is insensitive to variations Hubble constant, Hq, modulo an 
overall shift due to the change in the horizon size, rjQ. The WMAP likelihood of the corresponding 
Ci are ln£ = —477.79 and —478.83 at ho = 0.68 and 0.75, respectively. The variation of Baryon 
density around its middle value — 0.044 yields almost identical primordial power spectra with 
a minor change in likelihood -\n£ = -477.97 and -478.53 for 17b = 0.040 and = 0.048, 
respectively. The variation of JIa (or equivalently, flm, for fiat models) modulo a shift due to 
change in rjo, affects only the amplitude of the infra-red cut-off. The WMAP Hkelihood of the 
corresponding C/ are ln£ = —478.93 and —478.25 for fl^ = 0.69 and SIa = 0.77, respectively. 

We note that the likelihood at the central values are not necessarily the best. This suggests that, 
as discussed in the introduction, one should explore the entire space of cosmological parameters 
and compute the Hkelihood after optimizing with 'best' recovered primordial spectrum. Such an 
analysis is certainly possible, but has large demands on computational resources. Hence, we defer 
it to a future publication. 

Finally, we also consider the effect of a cosmological model with significant optical depth to 
reionization, t — 0.17, such as suggested by the angular power spectrum of the temperature- 
polarization (TE) cross correlation from WMAP fill . It has been pointed out that the estimated t 
increases when the primordial spectrum has an infra-red cut-off. We compute the primordial power 
spectrum for base cosmological models with t = 0.10, 0.17, 0.25. The standard RL deconvolution 
performs poorly for these models, and, the improved RL method was crucial for these models. 
Fig. IHl compares the recovered primordial spectra for early reionization models which all show an 
infra-red cutoff at the horizon scale [5fl| . 



IV. THEORETICAL IMPLICATIONS OF THE RECOVERED SPECTRA 



The direct recovery of the primordial power spectrum has revealed an infra-red cutoff of a very 
specific form. Model spectra with monotonic infra-red cutoff such as that in eq. H3.2|l do not 
improve the WMAP likeHhood significantly. While, to match the low value of the quadrupole, a 
very sharp cut-off (such as a 10, see figEJ is required, such a steep monotonic cut-off tends to 
pull down the power in the next few higher multipoles above the quadrupole and octopole as well. 
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FIG. 7: The panel of figures shows the robustness of the recovered P{k) for variations in the cosmological 
parameters. Each parameter is varied within the 1 a range indicated by the WMAP parameter estimates 
In Fig. (a), the Hubble constant ho = 0.68, ho = 0.71 and ho = 0.75 in the three curves. In Fig (b) the 
values of vacuum density — 0.69, = 0.73 and ~ 0.77 in the three curves. In fig. (c), the 
Baryonic density fis = 0.040, = 0.044 and = 0.048 in the three curves. The fig. (d) combines all 
the distinct curves in other figures to give a consolidated perspective on the dependence of the recovered 
spectrum on cosmological parameters. Note that the a;-axis is the wavenumber is scaled in units of the 
kfi = 2n/rio which reduces the scatter considerably in the curves for variations in Ho and D,a- 



Our recovered spectrum has a compensating excess which allows a steep cut-off to match the low 
quadrupole and octopole without suppressing the higher multipoles. Naively, one would think that 
a designer infra-red cutoff would 'cost' in the language of Bayesian evidence due to the introduction 
of extra parameters. That is is not necessarily so. An infra-red cut-off of the form we recover does 
not necessarily have more parameters than an infra-red cut-off of the form in eq. Ij^l2ll . Moreover, 
it is striking that the location of the cut-off is close to a well known scale - the horizon scale. 

In this section we show that infra-red cut-off of the form we recover arises from very simple 
scenarios in inflation. We explicitly mention two of them. Starobinsky 0| has shown that a kink 
(sharp change in the slope) in the inflaton potential can modulate the underlying primordial power 
spectrum Po{k) with a step like feature at a wavenumber kc 



P{k) 



Po{k)D{k,k„r 



1-"= [1 - 3(r - l)i((l - \) sin22/ + - cos2y) 



9 



1 



-(r - 1)S((1 + -)cos2y - - sin2y)] 
2 y 



(4.1) 
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FIG. 8: The recovered P{k) for different values of optical depth (r = 0.00, 0.10, 0.17 and 0.25). The width 
of smoothing used is different for different cases. This is the main cause of the small shift in the location 
of the infra-red cutoff. 



where we assume a power law Po{k) — As /c^^"% y — k/kc and r is the ratio of the slope dV/d(j) 
before and after kink in the infiaton potential. The power spectrum P{k) in eq. 114.111 has a step 
up (going to larger k ) for r < 1 and a step down feature for r > 1. An infra-red cut off is 
created with r < 1. Fig. El shows a spectrum with a Starobinsky step (eo. I4.1|l that can not only 
mimic the sharp infra-red cutoff but also produces the required bump after it. Tabled shows that 
a introducing an appropriate Starobinsky step gives a very good WMAP likehhood compared to 
eq. 113.211 . Besides the location of the break, kc, the Starobinsky step spectrum has only another 
parameter r that fixes the slope and the depth of the cut-off as well as the size of the bump. We 
have not systematically searched through the {r, kc} parameter space to arrive at a 'best-fit' model. 
Hence, it may be possible to get even better match to the WMAP data with Starobinsky breaks. 
Similar scenarios have been studied earlier [Soil and has also been pointed to in the post- WMAP 
literature J^S- Multiple scalar field infiation provide ample scope for generating features in the 
primordial spectrum Q and has been been invoked to model a sharp cut-off at horizon scale (see 
eg-, More exotic origin of an infra-red cut off in the scalar spectrum have also been 

investigated J53, 34, 35, 36, 32]. 

Another compelling theoretical scenario for generating a feature of the form we have recovered is 
well-known. It is well known that radiation, or matter dominated era prior to inflation does affect 
the primordial power spectrum on scales that 'exit the horizon' soon after the onset of inflation. 
For a pre-inflationary radiation dominated epoch the power spectrum was given by Vilenkin and 
Ford (VF) 



P{k) = As A;!-"^ — \e'^'y{l + 2iy) - 1 - 2y^\^ (4.2) 

where y — k/kc. Fig. El shows that the VF spectrum (ea. I4.2|l can also provide an infra-red cutoff 
with required bump after it. Tabled shows that a VF spectrum can give better WMAP hkelihood 
compared to model spectra of the form eq. 113. 2|l . The infra-red cut-off (oc k"^) here is not very 
sharp. However, if the epoch prior to inflation is dominated by matter with some other equation 
of state, the slope would be different. A more complete analysis may give rise to spectra closer to 
the kind we have recovered. The scale kc is set by the Hubble parameter at the onset of inflation. 
For this scenario to apphcable to our results, the fc/i-mode corresponding to the horizon scale must 
have crossed the Hubble radius very close to the onset of inflation. In a single scalar fleld inflation 
this would happen naturally (3^ . 
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FIG. 9: Comparison of our recovered P{k) (solid) with the predictions two simple theoretical scenarios 
that remarkably match the gross features of the infra-red cutoff in the recovered spectrum. The 'staro' 
curve is the primordial spectrum when the inflaton potential has a kink- a sharp, but rounded, change in 
slope U. Fine tuning is involved in locating the kink appropriately. The 'VF' curves are the modification 
to the power spectrum from a pre- inflationary (here, radiation dominated) epoch 3- This requires that the 
horizon scale, k^, exits the Hubble radius very soon after the onset of inflation. Although, it appears flne 
tuned there is corroborating support for this within single scalar fleld driven inflation [3g|- The theoretical 
P{k) leads to Ci that enhanced WMAP likelihood given in Table |l| The values of the parameters for the 
theoretical curves are given in the same table. 

Another possibility is that the infra-red cutoff arises due to non-trivial topology of the universe. 
A dodecahedral universe model that matches the low multipoles of WMAP angular spectrum has 
been proposed . In future work we would Hke to check whether the recovery of a discrete initial 
spectrum with our method appears similar to the spectrum we recovered here. However, non-trivial 
cosmic topology is expected to also violate the statistical isotropy of the CMB anisotropy and give 
rise of correlation features which are potentially detectable \4(\ . 

V. DISCUSSION AND CONCLUSIONS 

The CMB anisotropy is usually expected to be statistically isotropic and Gaussian ji^l. In that 
case, the angular power spectrum of CMB anisotropy encodes all the information that may be ob- 
tained from the primary CMB anisotropy, in particular, the estimation of cosmological parameters. 
It is very important to note that cosmological parameters estimated from CMB anisotropy (and 
other similar observations of the perturbed universe) usually assume a simple parametric form of 
spectrum of primordial perturbations. It is clear, however, that estimation of cosmological parame- 
ters depends on the extent and nature of parameterization of the primordial (initial) perturbations 
included into the parameter space considered 

We proceed on a complimentary path of determining the primordial power spectrum directly 
from the CMB anisotropy for a set of cosmological parameters. Assuming the best fit cosmological 
parameters from WMAP, our method applied to the angular power spectrum measured by WMAP 
yields interesting deviations from scale invariance in the recovered primordial power spectrum. 
The recovered spectrum shows an infrared cut-off that is robust to small changes in the assumed 
cosmological parameters. The recovered spectrum points to the form of infra-red cut off that 
matches the low multipoles of WMAP. We also show that the such forms of infra-red cut-off can 
arise from simple well-known effects with inflation. It is important to recall that the angular power 
spectrum from the 'full' sky CMB anisotropy measurement by COBE-DMR also indicated an 
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TABLE I: The effective chi-square, Xeff = — 21n£, of the Ci corresponding our recovered spectrum is 
compared with a number of model primordial spectrum (with or without the infra-red cutoff. Limited 
attempt has been made to search for the best parameter values and the Xeff for the model spectra should 
be treated as indicative and are strictly upper bounds. 



Power spectrum 


Xeff = -^ml. 


h / 1- ^ 






= 4.(3 X iU Mpc. j 


Direct Recovered 


f\t^£l 'TCI 

yOD. /o 


U. / 1 


Flat Harrison Zeldovich 


1007.28 




Power Law 


978.60 




(Hs = 0.95) 






Exponential cutoff 


978.08 


0.64 


{us = 0.95, a = 3.35) 






Exponential cutoff 


977.84 


0.64 


{us = 0.95, a = 10) 






Starobinsky break 


973.86 


0.32 


[us = 0.95, r = 0.01 ) 






Vilenkin & Ford (VF-1) 


976.88 


0.43 


{n, = 0.95 ) 






Vilenkin & Ford (VF-II) 


978.66 


0.96 


{us = 0.95 ) 







"Interpretation of fcc depends on the form of model power spectrum. 
For the recovered spectrum a simple tangent hyperbolic fit was used. 

infra-red cutoff . Although we mostly emphasize the infrared cutoff, the final recovered spectrum 
shows a damped oscillatory feature after the infra-red break ('ringing'). It has been pointed out 
recentl y th at such oscillations improve x^ff and could be possibly a signature of trans-Planckian 
effects [13,EJ|- We have not assessed the significance and robustness of the features at large kjkh 
that may be consistent with features deduced in the analysis of recent redshift surveys [45] . 

Here we have considered a 'best fit' cosmology and some finite number of variations around 
it However, with large but reasonable computational resource, it is possible to explore a 

multi-dimensional space of cosmological parameters to compute a likeHhood at each point for an 
optimal (recovered) primordial spectrum. 

In this work, we have used only the angular power spectrum of WMAP measured temperature 
fiuctuations (TT). Once the E-polarization power spectrum (EE) is announced by the WMAP 
team, our method can be extended to include TE and EE angular power spectra in obtaining the 
primordial power spectrum. It is also possible extend our recovery of the primordial spectrum to 
larger wave-numbers by using the matter density power spectrum measured by large scale redshift 
surveys such as the Sloan Digital sky survey SDSS and 2degree Field surve y, m easurement from 
Ly-a absorption, and possibly, weak gravitational lensing in the near future p9ll5(il|. 

In the absence of a definitive and precise scenario for the generation of primordial perturbations, 
a direct inversion of the primordial spectrum is an extremely appealing approach made possible by 
the remarkable quality of the recent cosmological data, in particular, the anisotropy in the cosmic 
microwave background temperature. 
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APPENDIX A: SOME ASPECTS OF THE METHOD. 

In this appendix, we provide support and justification for some steps in our method for recovering 
the primordial power spectrum. In section lA II we demonstrate the robustness of the iterative 
Richardson-Lucy deconvolution to changes in the initial guess. The next section IIA 2|l discusses 
the advantage of the improved Richardson-Lucy method used in our work. In the last section ll A 3|l . 
we model the broad features of reference spectrum analytically and show that it is well understood. 

1. Effect of the initial guess in the Richardson-Lucy method 
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FIG. 10: The figure shows the raw power spectrum recovered using the Richardson-Lucy algorithm starting 
from three different initial guesses. The effect of the initial guess is negligible in the other region k space. 
As shown in the next section lA 31 the artefacts at low k and high k which are removed by the reference 
spectrum have a known dependence on the initial guess. 

The effect of the initial guess is negligible in the Richardson-Lucy method of deconvolution for our 
problem. We find that the result for various different initial guesses come very close to each other 
after a few iterations. Here we demonstrate the robustness of our method by deconvolving using 
the Ci from a test spectrum (with a step) using different initial guesses. In Fig. ITTIl deconvolved 
raw P{k) arising from different forms of the initial guess spectra 

P{k) = constant, P{k) = 1/k, P{k) = (Al) 

are shown. The effect of the initial guess is absent in the portion of the k space probed well by 
the kernel. The deconvolved spectra are almost identical for all these different initial guesses in 
that region of k space. As discussed in section IXsl the dependence of deconvolved spectrum at the 
small and large k ends is well understood and get largely removed when divided by the reference 
spectrum. In a previous work we have checked a large variety of initial guesses and concluded the 
RL method applied to our problem of recovering the primordial spectrum from the CMB anisotropy 
is independent of the initial guess. 

2. The Improved Richardson Lucy method 

In this section we show the advantage of the improved Richardson Lucy algorithm given by 
eq. Ij2.14|l over the standard method given by eq. Ij2.13|l . We compare the raw power spectrum 
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FIG. 11: The left panel plots the variation of of C;'*' obtained after ith iteration (w.r.t the binned 
WMAP spectrum, Cf) with increasing iterations for the Richardson-Lucy (RL) method and improved 
version (IRL) we present in this work. The panel on the right, plots the variation of XeS = — 21n£ of the 
C;**' given by the WMAP likelihood, £. In contrast to RL, in the IRL method Xcft converges with iteration 
and to a significantly lower value. 

P^'^\k) and the corresponding angular spectrum C^^^ obtained from the WMAP binned angular 
power spectrum, Cf using RL and IRL algorithm as the iterations progress. In Fig. the left 

panel shows the of c|*-* w.r.t Cf. As expected, the RL method leads to a lower (better) than 
the IRL method. However, we are interested in a recovering a P{k) which give C/ that has better 
Hkelihood with respect to WMAP data. The right panel of Fig. [Til plots the WMAP Hkelihood 
given in terms of xls — — 21njC. This figure justifies the 'improved' label to the IRL method. For 
the RL case, the value of XcS ^^en to bounce and increase at large iterations. (This bounce 
is more pronounced and happens at a lower iteration for certain cosmological model, e.g., one 
with high optical depth to reionization, r.) This is a reflection of the problem of the RL method 
fitting the noise in the data. In contrast, in the IRL case the xls converges with iterations and to 
significantly better (lower) value that for the RL case. 

3. Reference spectrum 

In this subsection, we show that the reference spectrum used to remove numerical artifacts 
from the raw power spectrum recovered by the deconvolution is analytically well understood. The 
reference spectrum reflects the sensitivity of the kernel |A(Z,fc)|^ to the k space. (Here, the over- 
bar in |A(Z, k)\ alludes to the binned version related to G{1, k) following eq. II2.7|| V In the regions 
of k space where the kernel probes P(k) weakly, there is scope for changing the power spectrum 
without changing the C;. In such degenerate regions, the Richardson Lucy method (both RL and 
IRL) pushes the power spectrum up to the extent possible changing the Ci . 

Fig. El shows that the strong features at the low k and high k are well fitted by analytically ob- 
tained power law forms. At the low wavenumbers {k/kh ^ 1), the CMB anisotropy is dominated by 
the Sachs- Wolfe effect. The sensitivity of the CMB anisotropy to the power at very small wavenum- 
bers is dominated by the lowest multipole as |A(Z,fc)|^ oc fc^'. This is the well known Grishchuk- 
Zeldovich effect When the quadrupole is included, \A{1, fc)|^ cx fc^ is the strongest probe of the 
P{k) at low wavenumbers. Hence, there is no change to Ci (here, chiefly the quadrupole) caused 
from this region of k space from the initial guess P'"' (fc) onward if P{k) cx P(0)(fc) for e > 0. 

In our work we have used an initial guess of the form p(°)(fc) cx 1/k^. Hence, the slope of P{k) at 
low k is driven to the ~ k~^ form shown in Fig. El 

At the large k end, the |A(Z, fc)p is governed by the power slope of the tail in the sensitivity of 
the highest few multipole bins shown plotted in Fig. El The slope of these tails are well fltted by 
a power law form /c^^-^. Using the same argument, it is clear that there will be no change to the 
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FIG. 12: The figure shows the Praw(fc) and Pici{k) obtained from the WMAP binned data Cf for an 
initial guess P^''\k) oc l/k^. The dashed straight line labeled corresponds to the analytical power law 
form {k~'^P^''\k) = k~^) that matches the identical fall in Praw(fc) and Prcf(fc) at low k. The dashed 
line matching the rise in Praw(fc) and Prcf(fc) at large k corresponds to the analytical power law form 
(fc^-2p('')(fc) = fc5-2) expected from the roughly k ^'^ tail of the kernels for the last few multipole bins. 



Ci caused from this region of k space from the initial guess onward if P{k) oc p(o)(fc). For 

the initial guess P*^°)(/c) oc the rise at high k end is driven to a k^'^ form. 

The division of reference spectrum used here is a simple numerical recipe for estimating and 
removing the artefacts of the deconvolution method for a given kernel. We have shown that the 
gross features of the Prof (fc) can be understood analytically. The finer features of Prcf (fc) are likely 
Hnked to the details of the structure of |A(/, fe)p and will be explored in future work. 
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